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Abstract 

We derived the formulae of central differentiation for the finding of the 
first and second derivatives of functions given in discrete points, with the 
number of points being arbitrary. The obtained formulae for the deriva- 
tive calculation do not require direct construction of the interpolating 
polynomial. As an example of the use of the developed method we cal- 
culated the first derivative of the function having known analytical value 
of the derivative. The result was examined in the limiting case of infinite 
number of points. We studied the spectral characteristics of the weight 
coefficients sequence of the numerical differentiation formulae. The per- 
formed investigation enabled one to analyze the accuracy of the numerical 
differentiation carried out with the use of the developed technique. 

Mathematics Subject Classification: Primary 65D25; Secondary 65T50 

1 Introduction 

In solving many mathematical and physical problems by means of numerical 
methods one is often challenged to seek derivatives of various functions given in 
discrete points. In such cases, when it is difficult or impossible to take derivative 
of a function analytically one resorts to numerical differentiation. 

It should be noted that there exists a great deal of formulae and techniques 
of numerical differentiation (see, for instance, Ref. pQ). As a rule, the function 
in question f(x) is replaced with the easy-to-calculate function ip(x) and then it 
is approximately supposed that f'(x) ~ p'(x). The derivatives of higher orders 
are computed similarly. Therefore, in order to obtain numerical value of the 
derivative of the considered function it is necessary to indicate correctly the 
interpolating function ip(x). If the values of the function f(x) are known in s 
discrete points, the function ip(x) is usually taken as the polynomial of (s — l)th 
power. 

To find the derivative of functions having the intervals both quick and slow 
variation quasi-uniform nets are used (see Ref. |2])- This method has an ad- 



vantage since constant small mesh width is unfavorable in this case, because it 
leads to the strong enhancement of the function values table. 

The problem of the numerical differentiation accuracy is also of interest. The 
numerical differentiations formulae, taking into account the values of the con- 
sidered function both at x > xq and x < xq (xq is a point where the derivative 
is computed), are called central differentiation formulae. For instance, the for- 
mulae based on Stirling interpolating polynomial can be included in this class. 
Such formulae are known to have higher accuracy compared to the formulae, 
using unilateral values of a function 1 , i.e., for instance, at x > xq. 

The range of numerical differentiation formulae based on different interpo- 
lating polynomials is limited, as a rule, to finite points of interpolation. All 
available formulae known at the present moment are obtained for a certain con- 
crete limited number of interpolation points (see Refs. |H1E])- It can be explained 
by the fact that the procedure of the finding of the interpolating polynomial co- 
efficients in the case of the arbitrary number of interpolation points is quite 
awkward and requires formidable calculations. 

It is worth mentioning that the procedure of the numerical differentiation 
is incorrect. Indeed, in Ref. 2; it was shown that it is possible to select such 
decreasing error of the function in question which results in the unlimited growth 
of the error in its first derivative. 

Some recent publications devoted to the numerical differentiation problem 
should be mentioned (see, e.g., Ref. [5]). In this work the finite difference 
formulae for real functions on one dimensional grids with arbitrary spacing were 
considered. 

The formulae of central differentiation for the finding of the first and the 
second derivatives of the functions given in (2n + 1) discrete points are derived 
in this paper. The number of interpolation points is taken to be arbitrary. The 
obtained formulae for the derivatives calculation do not require direct construc- 
tion of the interpolating polynomial. As an example of the use of the developed 
method we calculate the first derivative of the function y(x) — sins. The ob- 
tained result is studied in the limiting case n — > oo. We examine the spectral 
characteristics of the weight coefficients sequence of the numerical differentiation 
formulae for the different number of the interpolation points. The performed 
analysis can be applied to the studying of the accuracy of the numerical differ- 
entiation technique developed in this work. It is found that the derived formulae 
of numerical differentiation have a high accuracy in a very wide range of spatial 
frequencies. 

The formulae using, for example, Newton interpolating polynomial are attributed to this 
class of numerical differentiation formulae. 



2 Formulae for approximate values of the first 
and the second derivatives 

Without the restriction of generality we suppose that the derivative is taken in 
the zero point, i.e. xq = 0. Let us consider the function f(x) given in equidistant 
points x m = ±mh, where m = 0, . . . , n and h is the constant value. We can 
pass the interpolating polynomial of the 2nth power through these points 



P2n(x) = ]T 



c k x 
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the values of the function in points of interpolation f m = f(x m ) coinciding with 
the values of the interpolating polynomial in these points: P2n{x m ) = fm- Let 
us define as d m the differences of the values of the function f(x) in diametrically 
opposite points x m and X- m , i.e. d m = f m — /_ m . We can present d m in the 
form 



(2.2) 
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To find the coefficients c^k+i, k — 0, ...to — 1, we have gotten the system of 
inhomogeneous linear equations with the given free terms d m . It will be shown 
below that this system has the single solution. 

We will seek the solution of the system [Eq. 12. 2|) ] in the following way 
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where a^ +1 ^ (to) are the undetermined coefficients satisfying the condition 
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Sik, l,k = 0,...,n — l. 



(2.4) 



Thus, the system of equations [Eq. Ij2.2|l ] is reduced to the equivalent, but more 
simple system [Eq. (|2.4|l ]. in which for each fixed number k = 0, . . . , n — 1 it is 
necessary to find the coefficients am l+1 \n). 

Let us resolve the system of equations [Eq. (|2.4I) ] according to the Cramer's 
rule: 
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In Eq. I|2.6fl we used the formula for the calculation of the Vandermonde de- 
terminant. From Eq. (|2.6I) it follows that the determinant of the system of 
equations [Eq. I|2.4fi ] is not equal to zero, i.e. the system of equations [Eq. 12.2(1 ] 
has the single solution. 

The most simple expression for Am +1 '(n) is obtained in the case of I = 
that corresponds to a calculation of the first-order derivative 

~* 2 )- (2-7) 
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From Eq. (|2.5|) as well as taking into account Eqs. (|2.6f) and (|2.7|) we get the 

expression for the coefficients otm (n) 

a ( i\n) = — \- (2.8) 

where 

n / 2 \ 

*m(n)= n (i-pj- (2-9) 

It should be noted that one can similarly get the expression for the coefficients 
otm 1 (n) which is presented in the following way 
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Taking into account Eqs. H2.lfl ~ H2.3fl we finally get the formula for the first 
derivative of the function f{x) 

1 ™ 

f'(0) « ^(0) = ^ E oWWifm - f-m). (2.10) 
rn—1 

The algorithm for the computation of the coefficients ajv i n ) is presented in the 
appendix and the results for the certain concrete number of the interpola- 
tion points (2n + 1) are given in Tab.^ Note that the expression for the first 
derivative obtained by this method coincides with the value presented in the 
Refs. [HI for n = 1,2 that corresponds to three and five points of interpola- 
tion. However, technique developed in this article allows one to calculate the 
coefficients a.m (n), and hence the first derivative, for any value of n. 
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-1/2772 



Table 1: Values of the coefficients a m (n). 



Similar formula can be obtained for the calculation of the second derivative. 
We give without proof corresponding expression 
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and the product TT m (n) is introduced in Eq. (|2.9() . 

As an example of the use of the obtained central differentiation formulae we 
will compute the first derivative of the function y(x) = sin x at x — 0. Let us 
set the value of the mesh width h equal to tv/2. Notice that, as a rule, the less 
the mesh width h the more exact result numerical differentiation gives. We have 
chosen rather big value of h. The Eq. I|2.10|l for this case takes the form 

2 ™ 

y'W = -J2(- i r4ll + i(n). (2.13) 

m=0 

In Eq. (|2.13() we take that d m = if m is an even number, and d m — ±2 if m 
is an odd number. 

Let us study the obtained result in the limiting case n — * oo. First, it is 
necessary to calculate the value of the product w m (n) within the limit n —> oo 
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Here we used the known value of infinite product 
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Using Eqs. I|2.8[) and Q2.14p . we find that the expression for the coefficients 
otm (n) within the limit n — > oo is represented in the following way 

a«= lim Q « (n) = ( _ ir+ il. (2 . 15) 

n— »oo 77j 

Now it is easy to complete the studying of Eq. I|2.13[) . Substituting the result 
from Eq. (|2.15|) to Eq. I|2.13[) and using the known value of infinite series we get 
that 

yKJ n ^ 2m + 1 

m— 

Thus, it is shown that the method of the derivatives finding, developed in this 
paper, gives for the function y(x) — sinx the value of the first derivative which 
coincides with the exact analytical one even at rather crude mesh width. 



3 Spectral characteristics of weight coefficients 
sequences 

In the section [21 of the present work we derived the formula for the finding of 
the first derivative of the function f(x) at x = 0. This result can be easily 
generalized for the case of the arbitrary point x — kh. If we set that a^] n (n) = 

—a£)(n), and moreover supposing that a™ (n) = for m = and m > n (see 
Tab.^), then in the considered case Eq. I|2.10|l reads as follows 

f(kh) « PL(kh) = ^ (3-1) 

m 

Here the summing is taken over all range of the function involved: f(kh). For 
instance, if the values of the function are set on the limited equidistant collection 
of elements N, then Eq. (|3 . II) can be rewritten in the form 

N-i 

771=0 

It is worth noticing that in Eq. (|3.2f> we used the periodicity condition of the 
weight coefficients 

Fig- presents the example of the weight coefficients of the differentiating se- 
quence Om(l). Thus the first derivative computation of the function f(x) at 
the points x — kh, k = 0, 1, . . . , N — 1, is reduced to the procedure of the cal- 
culation of the mutual correlation function between the finite sequences ot£) (n) 
and f m . 



taf ) (l) = l 



i 

i A r - 2 A' - 1 

-* • • t 



^i(l) = -l' 



Figure 1: The coefficients a™ (1). 



It is known (see, e.g., Ref. 6 ) that if a function satisfies the Dirichlet con- 
ditions in the interval (—1,1), then it can be expanded into the Fourier series 
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where the expansion coefficients are presented in the way 
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If the first derivative f'(x) satisfies the analogous conditions as the function 
f(x), then the following expression will be valid 



f( x )= X! |*t"} Cfc6xp (*~r a 
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Therefore, form Eqs. (|3.3(l and (|3.4(l it follows that the differentiation proce- 
dure is the linear filter with the frequency characteristic: &i(k) = ik{ir/l) 0. 
Similarly we receive for the second derivative 

/"OO = ]C j - (t) | Cfc6xp ( K i T x 



In this case the the frequency characteristic of the corresponding filter has the 
form: £ 2 (fc) = -k 2 {n/l) 2 . 

According to Wiener-Khinchin theorem (see, e.g., Ref. [7]) the mutual cor- 
relation function between the two finite sequences can be calculated with the 
help of the inverse Fourier transform of the mutual spectrum of the considered 
sequences. Thus, if we define that 



m=0 



3.5 



3 
2.5 
£ 2 

1 

0.5 



"0 200 400 600 800 1000 

r 

Figure 2: The spectra of various sequences a„ (n) at N — 2000. 
is the complex spectrum of the differentiating sequence a'm (n) , and 

c r=Yl /™ eX P ( 
m=0 ^ 

is the spectrum of the function f(x), then it follows from Eq. I|3.2[l that 

1 N ~ 1 / 2ir \ 

= 2h S Crft(r)expU—kr\ , (3.5) 

where /3*(r) is the complex conjugated quantity with respect to (3\{r). 

Comparing Eqs. I|3.4|l and 1|3.5[) we obtain that the accuracy of the numerical 
differentiation performed with the use of the various types of the sequences 
&rn {n) is characterized by the closeness of imaginary parts of their spectra to 
the linearly growing sequence yi(r) = 2nr/N. 

The spectra of the sequences a'm (ri) are depicted in Fig. [2] for the various 
values of n at N = 2000. It can be seen from this figure that for n = 1, i.e. 
for the sequence shown in Fig. ^ the imaginary part of the spectrum is the 
branch of the function sin(27rr/iV). The linearity condition is satisfied only in 
the vicinity of zero and N/2. However, at n = 11 the spectrum practically does 
not differ from the linear one up to r w N/2. The more close to linear one is 
the spectrum of the sequence am (21). 

The difference between the imaginary parts of the spectra of the sequences 
Qm (n) and the linearly growing sequence yi(r) = 2nr/N are presented in Fig.|5J 
The computations have been performed with the accuracy up to 10 -15 , thus the 
reliable results at n = 11 have been obtained for r > 150, and at n = 21 for 
r > 300. The presented results demonstrate the high accuracy of the numerical 
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Figure 3: The function tfi(r) = 5m(ft*(r)) — ?/i (r) versus r for different n. 

differentiation carried out with the help of the sequences am (n) in the wide 
range of the spatial frequencies. 

Now let us briefly consider the sequences for the calculation of the second 
derivative am (n), which are given in Eq. (|2.12(l . Their spectral properties 
can be obtained in the similar manner as we have done it for the case of the 
sequences am\n) and therefore we just present the final results. The spectra of 
the sequences am (n) are shown in Fig. 0] It follows from this figure that the 
closeness of the corresponding spectrum to the parabola 2/2O") = — (27rr/iV) 2 in 
the case of n = 1 exists only in the vicinity of zero. The spectra at n = 11 
and n = 21 are close to function y2(r) in a wider range of r (r < 750 and 
r < 800 respectively). The difference between the real parts of the spectra 
of the sequences am\n) and the parabola yz(r) = —{2-Kr/N) 2 are depicted in 
Fig.JSJin the logarithmic scale. This figure again demonstrates the high accuracy 
of the second derivative computation with the use of the sequences a m (n). 

4 Conclusion 

In conclusion we note that the method of central differentiation formulae finding 
has been developed in this article. The elaborated technique does not require 
direct construction of the interpolating polynomial. We have derived simple and 
convenient expressions for the first and the second derivatives [Eqs. (|2.10JI and 
H2.11f> ] of the function given in (2n + l) discrete points. The number n was taken 
to be arbitrary. In contrast to the results of the Ref. [5], where the recursion 
relations for the calculation of the weight coefficient being used in numerical 
differentiation formulae were considered, in the present work the expressions 
for the considered weight coefficients have been derived in the explicit form for 
the arbitrary number of interpolation points. As an example of the use of the 
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Figure 5: The function S 2 (r) = Ut(ffi(r)) - y 2 (r) versus r for different n. 



developed method we have calculated the first derivative of the function y(x) — 
sin x. The obtained result has been studied in the limiting case ti-*oo. We have 
examined the spectral characteristics of the weight coefficients sequence of the 
numerical differentiation formulae for the different number of the interpolation 
points. The performed analysis has allowed one to study the accuracy of the 
numerical differentiation carried out with the help of the developed method. It 
has been found that the derived formulae of numerical differentiation posses the 
high accuracy in a rather wide range of the spatial frequencies. As it has been 
shown in this paper, the formulae for the derivatives finding gave correct results 
in the case of large number of interpolation points. Thus, the developed method 
can be useful in lattice simulation of quantum fields 8 . To get the exact results 
at calculations on lattices one has to use nets with the big number of points. 
Derivatives which one encounters in theories of quantum fields, as a rule, do not 
exceed the second order. Therefore, the formulae obtained in this article could 
be of use in carrying out mentioned above research. 
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A Algorithm for computation of weight coeffi- 
cients am{n) 

In this appendix we present the algorithm for the computation of the coefficients 

&m\n) on the MATLAB 6.5 programming language. 

N = 2000; N is the size of the array 

n = 11; n is the order of the differentiating sequence 

a = zeros(l, N); a is the array of the weight coefficients 

fci = 2; k 2 = N; 
for 77i — 1 : n 

n = l; 

for k = 1 : n 

if k == 777 

T2 = l; 

else 

r 2 = 1 - (m/k) 2 ; 
end 

n = n * r-x\ 

end 

r\ = 1/(2 * ri * m); 

a(l, ki) = — n; k\ = k\ + 1; 

a(l, k 2 ) = n; k 2 = k 2 - 1; 

end 
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